function [vals, conds] = eval_range(sF, v, varargin)

% get parameters
dimensions = size(v);
v = v(:);
N = size(v, 1);
vals = zeros(N, numel(sF));
conds = zeros(N, 1);
 
% get the neighbors 
ind = sF.nodes.find(v, sF.delta); 
nn = sum(ind, 2);

% for points with too less neighbors, we instead choose the sF.dim nearest ones
% NOTE: we choose more neighbors than only the sF.dim nearest ones, since the
% expectation of the lebesgue constant is infinite in this setting
I = nn < sF.dim;
if (sum(I) > 0)
  warning(sprintf( ...
    ['Some centers did not have sufficiently many neighbors. \n' ...
    '\t In this case the numer of neighbors was set to the dimension of the ansatz space.']));
  
  % evaluate the critical nodes via knn-search instead of rangesearch
  nn_original = sF.nn;
  sF.nn = sF.dim;
  if (nargout == 2)
    [vals(I,:), conds(I)] = sF.eval(v.subSet(I));
  else
    vals(I,:) = sF.eval(v.subSet(I));
  end
  sF.nn = nn_original;

  if (sum(I) == N)
    return;
  end
end

% continue with the points that have sufficiently many neighbors
J = ~I;
v = v.subSet(J);
N = sum(J);
[ind, dist] = sF.nodes.find(v, sF.delta);
% if optimal subsampling is set to true, we can now fall back to the eval_knn case 
%   where all neighborhoods have the same size (the dim of the ansatz space) 
if (sF.subsample == true)
  ind = sF.find_optimal_subset(ind, v, varargin{:});
end

[grid_id, v_id] = find(ind');
nn = sum(ind, 2);

if (sF.subsample == true)
  dist = angle(v.subSet(v_id), sF.nodes.subSet(grid_id));
  dist = sparse(v_id, grid_id, dist, N, numel(sF.nodes));
end


% the index vector col_id helps to construct the (sF.dim x N) matrix G, which
% holds the values of the basis functions at all neighbors of all centers from v
% col_id skips entries, whenever a center has not nn_max many neighbors 
nn_total = sum(nn);
nn_max = max(nn);
start_id = cumsum(nn(1:N-1)) + 1;
temp = ones(nn_total, 1);
temp(start_id) = 1 - nn(1:N-1);
temp = cumsum(temp);
col_id = (v_id-1) * nn_max + temp;

% TODO: nn_max might be much larger than mean(nn) at very few occations
%   ==> compute in batches of similar nn for less ram usage

% compute for every center from v the matrix of all basis functions evaluated at
% all neighbors of this center 
G = zeros(sF.dim, nn_max * N);
% evaluate the basis functions on the nodes
if (~sF.centered)
  % choose faster way between computing all values and reusing them or
  % computing values on fibgrid(grid_id)
  if nn_total > numel(sF.nodes.x)
    basis_on_grid = eval_basis_functions(sF);
    G(:, col_id) = basis_on_grid(grid_id, :).';
  else
    G(:, col_id) = eval_basis_functions(sF, sF.nodes(grid_id)).';
  end

  % odd basis functions may clash with antipodal option, since (-v) = -p(v)
  % thus make sure to use the representer which is closer to the center
  if (mod(sF.degree, 2) > 0)
    I = sum(v.subSet(v_id).xyz .* sF.nodes.subSet(grid_id).xyz, 2) < 0;
    G(:,I) = G(:,I) * (-1);
  end

  basis_in_v = eval_basis_functions(sF, v);
else
  % compute the rotations that shift each element of v into the north pole
  rot = rotation.map(v, vector3d.Z);
  rot = rot(v_id);
  rotneighbors = rot .* sF.nodes(grid_id);

  % determine which basis to use and evaluate it on the grid and on v
  basis_on_grid = eval_basis_functions(sF, rotneighbors);
  basis_in_pole = eval_basis_functions(sF, vector3d.Z);
  
  basis_in_v = repmat(basis_in_pole, N, 1);
  G(:, col_id) = basis_on_grid.';
end

% dont solve the normal equations G'WGc = G'Wf (like cond(G)^2)
% rather let matlab directly find min norm solution of sqrt(W) * (Gc-f)
% internally this uses QR and we end up with only cond(G)

% compute the weights
weights = zeros(N * nn_max, 1);
% dist(find(ind)) instead of nonzeros(dist), since elements of v might be
%   contained in sF.nodes ==> distance 0
weights(col_id) = sF.w(dist(ind > 0) / sF.delta);

B = G .* sqrt(weights');
B_book = pagetranspose(reshape(B, sF.dim, nn_max, N)); 

% compute scaling factors (norms of columns of G_times_W_book)
S_book = sqrt(sum(abs(B_book).^2, 1));

% set up right hand side
f = zeros(N * nn_max, numel(sF));
f(col_id,:) = sF.values(grid_id,:);
fw_book = permute(reshape((sqrt(weights) .* f).', numel(sF), nn_max, N), [2 1 3]);

% compute the generating functions
c_book = pagemldivide(B_book ./ S_book, fw_book) ./ pagetranspose(S_book);
vals(J,:) = permute(sum(basis_in_v .* permute(c_book, [3 1 2]), 2), [1 3 2]);

if isalmostreal(sF.values)
  vals = real(vals); 
end

if nargout == 2
  eigsJ = pagesvd(B_book ./ S_book);
  condsJ = eigsJ(1,:,:) ./ eigsJ(sF.dim,:,:);
  conds(J) = condsJ(:);
  conds = reshape(conds, dimensions);
end

end
